# PROLOG   ##############################################################################

# PROJECT: Priority Effects Experiment 
# PURPOSE: Determine how species arrival via seed additions affects the community assembly of restored prairie plant communities
#
# AUTHOR:  Katherine Wynne (Wynnekat@msu.edu)
# COLLAB:  L. Sullivan
# CREATED: 06 December 2022
#
# FILES:   1) inner_cover_summer_2022.xlsx - vegetative cover taken in the Summer (end of June)
#          2) inner_cover_fall _2022.xlsx - vegetative cover taken at peak biomass (end ofAugust)
#          3) Checklist of Missouri Flora.xlsx - contains the information about all plants found in MO
#          4) Cleaned_PE_Detailed_Species_List_2022.xlsx - dataset containing the information about only                                                                plants found in the study
#         
#LAST UPDATED: 06 September 2023
# PROLOG   ##############################################################################

Treatment information & Setup

NON = None (no species seeded)

SIM = Simulataneous (36 species seeded on March 22nd, 2021)

LE = Lumped early (18 early-dispersing species seeded on March 23rd, 2021 then 18 late-dispersing species seeded on September 5th, 2021)

LL = Lumped late (18 late-dispersing species seeded on March 23rd, 2021 then 18 early-dispersing species seeded on September 5th, 2021)

NAT = Natural; species added according to peak dispersal date

    - May 24th 2021 (VIOPED, PACPLA)
    - June 4th 2021 (SISCAM, SPHOBT, CORLAN)
    - June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC)
    - July 11th 2021 (CORPAL, DODMEA)
    - July 25th 2021 (KOEMAC, AMOCAN)
    - August 2nd 2021 (LINSUL, MELVIR)
    - August 21st 2021 (DALCAN, DALPUR, ACHMIL)
    - September 17th 2021 (CROSAG)
    - October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR)
    - October 18th 2021 (BIDARI, HELMOL, LESCAP)
    - November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC)
    - November 19th 2021 (CORTRI, PYCTEN, SOLRIG)
  

Load libraries

Import Datasets

Dataset cleaning

2021

### Remove unnecessary columns

#### Get rid of the notes, unknown, and senensced columns

### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]


#### change log
 # - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct.  (Sep 6, 2023)
 # - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)


### Get rid of EMP

### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")

2022

2023

# Make a unique identifier for year and season

### 2021

#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))


### 2022

# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)

### 2023

#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))

#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)

inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)
# Join years together

full_2021_data <- inner_cover_fall_2021

full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)
## Joining with `by = join_by(Season, Year, Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name)`
full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
# Join all years together

### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)
## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`
# lump certain taxa together 

### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)

## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP

for(i in 1:nrow(full_year_data)) {
  if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "LEPSPP"){full_year_data [i,3] <- "LEPVIR"}
  if(full_year_data[i,3] == "ECHSPP"){full_year_data [i,3] <- "ECHCRU"}
  
  # Make all cover data that was less than 1 = to 1 (to make the data discrete)
  if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
  
}
inner_cover_max <- full_year_data  %>%
      group_by(Block, Treatment, Year, SPP6) %>%
      summarise(max_cover=max(Percent_Cover))
## `summarise()` has grouped output by 'Block', 'Treatment', 'Year'. You can
## override using the `.groups` argument.
##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice


### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")
### Remove Litter
inner_cover_max_only  <- subset(inner_cover_max_only, SPP6 != "LITTER")

head(inner_cover_max_only) 
## # A tibble: 6 × 5
## # Groups:   Block, Treatment, Year [1]
##   Block Treatment Year  SPP6   max_cover
##   <dbl> <chr>     <chr> <chr>      <dbl>
## 1     1 LE        2021  ACAVIR         1
## 2     1 LE        2021  ACHMIL         1
## 3     1 LE        2021  AMATUB         1
## 4     1 LE        2021  AMBART        20
## 5     1 LE        2021  BARVUL         1
## 6     1 LE        2021  BROJAP         1

Bare/Litter

### Bare dataset

inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")

### Litter dataset

inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")


### Manipulate 2021 dataset to match 2022 and 2023

#### Bare Ground

bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))

for(i in 1:nrow(bare_2021)) {
  if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}

#### Litter 

litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))


for(i in 1:nrow(litter_2021)) {
  if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}


### Full_join datasets

inner_bare_cover <- full_join(inner_bare_cover, bare_2021)
## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`
inner_litter_cover <- full_join(inner_litter_cover, litter_2021)
## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`
### Add species information to the dataset

#### N/A is N = native, A = non-native, G = Genus

inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)
## Joining with `by = join_by(SPP6)`
#Remove unnecessary columns
inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]

Exploratory Analysis of Diversity

Summary Table of Diversity Indices

Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.

Analyzing diversity

Mean C

Mean C - Model assumptions

Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. SIM had the most “even” distribution of C values and almost every value was represented in this treatment.

## Warning in geom_density(bins = 30, alpha = 0.4): Ignoring unknown parameters:
## `bins`

# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT

# only one extreme value in LE 2021 (likely because of DALCAN and VIOPED)
Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Mean_C)
## # A tibble: 5 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <dbl>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 NAT       2023      4               18    2.44   2.69 TRUE       FALSE     
## 2 NAT       2023      5               12    1.73   0.5  TRUE       FALSE     
## 3 NON       2023      6               19    1.99   1.18 TRUE       FALSE     
## 4 SIM       2022      5               20    1.83   3.13 TRUE       FALSE     
## 5 SIM       2023      2               27    2.41   3.42 TRUE       FALSE
# Mean C fits a normal distribution

Summary_table  %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Mean_C)
## # A tibble: 10 × 5
##    Treatment Year  variable statistic      p
##    <chr>     <chr> <chr>        <dbl>  <dbl>
##  1 LE        2022  Mean_C       0.926 0.548 
##  2 LE        2023  Mean_C       0.953 0.767 
##  3 LL        2022  Mean_C       0.890 0.317 
##  4 LL        2023  Mean_C       0.861 0.193 
##  5 NAT       2022  Mean_C       0.901 0.379 
##  6 NAT       2023  Mean_C       0.929 0.572 
##  7 NON       2022  Mean_C       0.928 0.566 
##  8 NON       2023  Mean_C       0.882 0.280 
##  9 SIM       2022  Mean_C       0.816 0.0809
## 10 SIM       2023  Mean_C       0.879 0.263
ggplot(data = Summary_table , aes(x= Mean_C)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table , "Mean_C", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

Mean C - Model fitting

Mean C values were significantly higher in all treatments compared to the NON treatment.

#Dropped the mixed model due to singular fit -> the variance of the random effect block is estimated close to zero and does not further inform the data. 

Summary_table$Block <- as.factor(Summary_table$Block)

# Normal interaction model

#Mean_C.mod.lmer.interaction <- lmer(Mean_C~Treatment*Year+(1|Block), data =Summary_table)
#AIC(Mean_C.mod.lmer.interaction)


# Normal - additive model (has the lowest AIC)


Mean_C.mod.lm.additive <- lm(Mean_C~Treatment+Year, data =Summary_table)

## Summary results
summary(Mean_C.mod.lm.additive)
## 
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07124 -0.32383 -0.07948  0.28266  1.12107 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.6602     0.1608  10.327 2.16e-14 ***
## TreatmentLL    0.5149     0.2075   2.481 0.016251 *  
## TreatmentNAT  -0.3549     0.2075  -1.710 0.093003 .  
## TreatmentNON  -1.3445     0.2075  -6.479 2.89e-08 ***
## TreatmentSIM   0.7983     0.2075   3.847 0.000318 ***
## Year2023       0.2660     0.1313   2.026 0.047691 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5084 on 54 degrees of freedom
## Multiple R-squared:  0.7134, Adjusted R-squared:  0.6869 
## F-statistic: 26.88 on 5 and 54 DF,  p-value: 1.544e-13
## AIC test
AIC(Mean_C.mod.lm.additive)
## [1] 96.76091
## ANOVA type 3
anova(Mean_C.mod.lm.additive)
## Analysis of Variance Table
## 
## Response: Mean_C
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Treatment  4 33.674  8.4186 32.5772 8.088e-14 ***
## Year       1  1.061  1.0610  4.1057   0.04769 *  
## Residuals 54 13.955  0.2584                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## I removed the random effect block because it had no effect on the model fit and explained no additional variance. (chose the most parsimonius model)

# r.squaredGLMM(Mean_C.mod.lmer.additive)


## Inclusion of the random effect marginally increased model fit 



est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.additive, "Treatment")


#Posthoc test

## Treatment
pairs(est.lm.treatment, adjust = "tukey")
##  contrast  estimate    SE df t.ratio p.value
##  LE - LL     -0.515 0.208 54  -2.481  0.1102
##  LE - NAT     0.355 0.208 54   1.710  0.4366
##  LE - NON     1.345 0.208 54   6.479  <.0001
##  LE - SIM    -0.798 0.208 54  -3.847  0.0028
##  LL - NAT     0.870 0.208 54   4.191  0.0009
##  LL - NON     1.859 0.208 54   8.960  <.0001
##  LL - SIM    -0.283 0.208 54  -1.366  0.6519
##  NAT - NON    0.990 0.208 54   4.769  0.0001
##  NAT - SIM   -1.153 0.208 54  -5.557  <.0001
##  NON - SIM   -2.143 0.208 54 -10.325  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 5 estimates
# Residual plot looks decent
plot(Mean_C.mod.lm.additive)

Mean C - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.c <- Summary_table_mean$C_Mean - Summary_table_mean$C_SD
upper.c <- Summary_table_mean$C_Mean + Summary_table_mean$C_SD

pd <- position_dodge(width = 0.8)

C.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = C_Mean, color = Treatment))+
  geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.c , ymax = upper.c), width = 0.5, position = pd) +
  labs(y = "Mean C", x = "Year" )+
  theme_classic() +
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))

C.plot

Species Richness

SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.

SR - Model assumptions
# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT

Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Species_richness)
## # A tibble: 5 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <fct>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 LE        2023  3                   15    2.10   2.2  TRUE       FALSE     
## 2 LL        2022  2                   25    2.31   2.06 TRUE       FALSE     
## 3 NAT       2022  4                   22    2.18   2.31 TRUE       TRUE      
## 4 NAT       2022  6                   15    2.10   1.09 TRUE       FALSE     
## 5 NON       2023  6                   19    1.99   1.18 TRUE       FALSE
# Species richness fits a normal distribution

Summary_table %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Species_richness)
## # A tibble: 10 × 5
##    Treatment Year  variable         statistic      p
##    <chr>     <chr> <chr>                <dbl>  <dbl>
##  1 LE        2022  Species_richness     0.996 0.998 
##  2 LE        2023  Species_richness     0.847 0.149 
##  3 LL        2022  Species_richness     0.814 0.0778
##  4 LL        2023  Species_richness     0.825 0.0969
##  5 NAT       2022  Species_richness     0.836 0.122 
##  6 NAT       2023  Species_richness     0.884 0.286 
##  7 NON       2022  Species_richness     0.960 0.817 
##  8 NON       2023  Species_richness     0.955 0.784 
##  9 SIM       2022  Species_richness     0.893 0.332 
## 10 SIM       2023  Species_richness     0.873 0.238
ggplot(data = Summary_table, aes(x= Species_richness)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table, "Species_richness", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

SR - Model fitting
# Model Fitting

## Despite being count data, the distribution of Species richness is remarkabley balanced and not skewed. Model comparison using AIC reveals that the Normal model works better than the Poisson model. Model results also differ likely because of under dispersion. The more complicated Quasipoisson model that accounts for under dispersion produces similar results to the Normal mixed model. Therefore, I'm going to go with the less complicated Normal model.


SR.mod.normal <- lmer(Species_richness~Treatment+Year+(1|Block), data =Summary_table)
summary(SR.mod.normal)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 283.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3018 -0.7002  0.0875  0.6078  1.6946 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 2.378    1.542   
##  Residual             7.355    2.712   
## Number of obs: 60, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   19.9333     1.0639 20.8212  18.736 1.64e-14 ***
## TreatmentLL    0.3333     1.1072 49.0000   0.301 0.764634    
## TreatmentNAT  -1.7500     1.1072 49.0000  -1.581 0.120398    
## TreatmentNON  -4.3333     1.1072 49.0000  -3.914 0.000280 ***
## TreatmentSIM   4.2500     1.1072 49.0000   3.839 0.000355 ***
## Year2023      -1.7000     0.7002 49.0000  -2.428 0.018911 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.520                            
## TreatmntNAT -0.520  0.500                     
## TreatmntNON -0.520  0.500  0.500              
## TreatmntSIM -0.520  0.500  0.500  0.500       
## Year2023    -0.329  0.000  0.000  0.000  0.000
## Treatment and Year significant
## Interaction between Treatment and Year not significant

anova(SR.mod.normal, type = 3)
## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment 474.77  118.69     4    49 16.1381 1.704e-08 ***
## Year       43.35   43.35     1    49  5.8941   0.01891 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Marginal R^2 (variance explained by ONLY fixed effects) = 0.4547589
## Conditional R^2 (variance explained by the entire model) = 0.5554578

## Inclusion of the random effect increased model fit 


r.squaredGLMM(SR.mod.normal)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.4743091 0.6027586
#Posthoc test

est.lmer.SR.treatment <- emmeans::emmeans(SR.mod.normal, "Treatment")


## Treatment
pairs(est.lmer.SR.treatment , adjust = "tukey")
##  contrast  estimate   SE df t.ratio p.value
##  LE - LL     -0.333 1.11 49  -0.301  0.9981
##  LE - NAT     1.750 1.11 49   1.581  0.5166
##  LE - NON     4.333 1.11 49   3.914  0.0025
##  LE - SIM    -4.250 1.11 49  -3.839  0.0031
##  LL - NAT     2.083 1.11 49   1.882  0.3406
##  LL - NON     4.667 1.11 49   4.215  0.0010
##  LL - SIM    -3.917 1.11 49  -3.538  0.0076
##  NAT - NON    2.583 1.11 49   2.333  0.1519
##  NAT - SIM   -6.000 1.11 49  -5.419  <.0001
##  NON - SIM   -8.583 1.11 49  -7.753  <.0001
## 
## Results are averaged over the levels of: Year 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates
# Poisson Model

# SR.mod.pois <- glmer(Species_richness~Treatment+Year+(1|Block),family = "poisson", data =Summary_table)
# summary(SR.mod.pois)
# Anova(SR.mod.pois, type = 3)

# Check dispersion

## Underdispersed! 0.659 
## Interaction between treatment and year are not significant

# Dispersion_glmer(SR.mod.pois)


# Not accounting for the underdispersion mattered a lot! .

## Used a quassipoisson model instead and year is significant

#m6 <-glmmPQL(Species_richness~Treatment+Year,
             # random = ~ 1 | Block,
             # family = quasipoisson(link='log'), 
             # data = Summary_table)

# summary(m6)
# Anova(m6, type = 3)
# Residual plot looks decent except for one large outlier
plot(SR.mod.normal)

SR - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.SR.c <- Summary_table_mean$SR_Mean - Summary_table_mean$SR_SD
upper.SR.c <- Summary_table_mean$SR_Mean + Summary_table_mean$SR_SD

pd <- position_dodge(width = 0.8)

SR.plot1 <- ggplot(data = Summary_table_mean, aes(x = Year, y = SR_Mean, color = Treatment))+
  geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.SR.c , ymax = upper.SR.c), width = 0.5, position = pd) +
  labs(y = "Total Species Richness", x = "Year" )+
  theme_classic() +
  ylim(10,30)+
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
    theme(legend.position="bottom")

SR.plot <- SR.plot1 +
    theme(legend.position="none")


leg <- get_legend(SR.plot1)

Shannon Diversity Index

SDI - Model assumptions

SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.

# Two extreme outliers (one in NAT 2021 and another in SIM 2021)
Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Shannon)
## # A tibble: 1 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <fct>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 NON       2022  3                   11    1.55  0.429 TRUE       FALSE
# Except for one point, Shannon diversity fits a normal distribution

Summary_table %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Shannon)
## # A tibble: 10 × 5
##    Treatment Year  variable statistic     p
##    <chr>     <chr> <chr>        <dbl> <dbl>
##  1 LE        2022  Shannon      0.971 0.901
##  2 LE        2023  Shannon      0.967 0.873
##  3 LL        2022  Shannon      0.977 0.936
##  4 LL        2023  Shannon      0.854 0.169
##  5 NAT       2022  Shannon      0.950 0.743
##  6 NAT       2023  Shannon      0.931 0.586
##  7 NON       2022  Shannon      0.848 0.151
##  8 NON       2023  Shannon      0.950 0.739
##  9 SIM       2022  Shannon      0.872 0.235
## 10 SIM       2023  Shannon      0.940 0.663
ggplot(data = Summary_table, aes(x= Shannon)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table, "Shannon", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

SDI - Model fitting
Shannon.lmer.mod <- lmer(Shannon~Treatment+Year+(1|Block), data =Summary_table)
summary(Shannon.lmer.mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 2.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60837 -0.47607  0.08046  0.51127  2.04855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 0.007995 0.08942 
##  Residual             0.042202 0.20543 
## Number of obs: 60, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   2.16123    0.07452 28.51199  29.003  < 2e-16 ***
## TreatmentLL   0.05159    0.08387 49.00000   0.615  0.54133    
## TreatmentNAT -0.12219    0.08387 49.00000  -1.457  0.15152    
## TreatmentNON -0.34644    0.08387 49.00000  -4.131  0.00014 ***
## TreatmentSIM  0.12917    0.08387 49.00000   1.540  0.12994    
## Year2023     -0.01536    0.05304 49.00000  -0.290  0.77333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.563                            
## TreatmntNAT -0.563  0.500                     
## TreatmntNON -0.563  0.500  0.500              
## TreatmntSIM -0.563  0.500  0.500  0.500       
## Year2023    -0.356  0.000  0.000  0.000  0.000
anova(Shannon.lmer.mod, type = 3)
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment 1.65266 0.41316     4    49  9.7902 6.691e-06 ***
## Year      0.00354 0.00354     1    49  0.0839    0.7733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# No significant interaction between treatment and year



## Marginal R^2 (variance explained by ONLY fixed effects) = 0.330704
## Conditional R^2 (variance explained by the entire model) = 0.3844412

## Inclusion of the random effect increased model fit 

r.squaredGLMM(Shannon.lmer.mod)
##           R2m       R2c
## [1,] 0.358652 0.4608049
# Post-Hoc comparison

pairs(emmeans::emmeans(Shannon.lmer.mod, "Treatment"))
##  contrast  estimate     SE df t.ratio p.value
##  LE - LL    -0.0516 0.0839 49  -0.615  0.9720
##  LE - NAT    0.1222 0.0839 49   1.457  0.5947
##  LE - NON    0.3464 0.0839 49   4.131  0.0013
##  LE - SIM   -0.1292 0.0839 49  -1.540  0.5420
##  LL - NAT    0.1738 0.0839 49   2.072  0.2486
##  LL - NON    0.3980 0.0839 49   4.746  0.0002
##  LL - SIM   -0.0776 0.0839 49  -0.925  0.8858
##  NAT - NON   0.2242 0.0839 49   2.674  0.0727
##  NAT - SIM  -0.2514 0.0839 49  -2.997  0.0331
##  NON - SIM  -0.4756 0.0839 49  -5.671  <.0001
## 
## Results are averaged over the levels of: Year 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates
# A little heteroscedastic (bowtie shaped but probably ok)
plot(Shannon.lmer.mod)

SDI - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.SDI.c <- Summary_table_mean$SDI_Mean - Summary_table_mean$SDI_SD
upper.SDI.c <- Summary_table_mean$SDI_Mean + Summary_table_mean$SDI_SD

pd <- position_dodge(width = 0.8)

SDI.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = SDI_Mean, color = Treatment))+
  geom_point(data =Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.SDI.c , ymax = upper.SDI.c), width = 0.5, position = pd) +
  labs(y = "SDI", x = "Year" )+
  theme_classic() +
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))

SDI.plot 

Plot Panel

Comparing seeded cover (in the inner 1 m^2 sampling area) for each treatment

note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)

Total Seeded Cover

TSC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year, 
##     family = binomial, data = Total_Cover_Tot)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.4760  -2.0224  -0.4559   1.5998   4.5960  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -2.11998    0.08881 -23.870  < 2e-16 ***
## TreatmentLL            0.94324    0.10412   9.059  < 2e-16 ***
## TreatmentNAT          -0.38112    0.13799  -2.762  0.00575 ** 
## TreatmentSIM           1.06270    0.10397  10.221  < 2e-16 ***
## Year2023               0.62783    0.11282   5.565 2.63e-08 ***
## TreatmentLL:Year2023  -0.41302    0.13685  -3.018  0.00254 ** 
## TreatmentNAT:Year2023 -0.36533    0.18174  -2.010  0.04441 *  
## TreatmentSIM:Year2023 -0.37777    0.13662  -2.765  0.00569 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 837.79  on 47  degrees of freedom
## Residual deviance: 293.23  on 40  degrees of freedom
## AIC: 557.58
## 
## Number of Fisher Scoring iterations: 5
## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                LR Chisq Df Pr(>Chisq)    
## Treatment       272.445  3  < 2.2e-16 ***
## Year             32.035  1  1.514e-08 ***
## Treatment:Year   10.305  3    0.01614 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast                    odds.ratio     SE  df null z.ratio p.value
##  LE Year2022 / LL Year2022        0.389 0.0405 Inf    1  -9.059  <.0001
##  LE Year2022 / NAT Year2022       1.464 0.2020 Inf    1   2.762  0.1048
##  LE Year2022 / SIM Year2022       0.346 0.0359 Inf    1 -10.221  <.0001
##  LE Year2022 / LE Year2023        0.534 0.0602 Inf    1  -5.565  <.0001
##  LE Year2022 / LL Year2023        0.314 0.0328 Inf    1 -11.076  <.0001
##  LE Year2022 / NAT Year2023       1.126 0.1469 Inf    1   0.909  0.9853
##  LE Year2022 / SIM Year2023       0.269 0.0281 Inf    1 -12.574  <.0001
##  LL Year2022 / NAT Year2022       3.760 0.4466 Inf    1  11.150  <.0001
##  LL Year2022 / SIM Year2022       0.887 0.0680 Inf    1  -1.558  0.7751
##  LL Year2022 / LE Year2023        1.371 0.1210 Inf    1   3.572  0.0085
##  LL Year2022 / LL Year2023        0.807 0.0625 Inf    1  -2.774  0.1016
##  LL Year2022 / NAT Year2023       2.892 0.3181 Inf    1   9.654  <.0001
##  LL Year2022 / SIM Year2023       0.691 0.0534 Inf    1  -4.784  <.0001
##  NAT Year2022 / SIM Year2022      0.236 0.0280 Inf    1 -12.169  <.0001
##  NAT Year2022 / LE Year2023       0.365 0.0461 Inf    1  -7.977  <.0001
##  NAT Year2022 / LL Year2023       0.215 0.0256 Inf    1 -12.917  <.0001
##  NAT Year2022 / NAT Year2023      0.769 0.1096 Inf    1  -1.842  0.5909
##  NAT Year2022 / SIM Year2023      0.184 0.0219 Inf    1 -14.231  <.0001
##  SIM Year2022 / LE Year2023       1.545 0.1361 Inf    1   4.935  <.0001
##  SIM Year2022 / LL Year2023       0.909 0.0702 Inf    1  -1.234  0.9218
##  SIM Year2022 / NAT Year2023      3.259 0.3580 Inf    1  10.753  <.0001
##  SIM Year2022 / SIM Year2023      0.779 0.0600 Inf    1  -3.246  0.0258
##  LE Year2023 / LL Year2023        0.588 0.0523 Inf    1  -5.971  <.0001
##  LE Year2023 / NAT Year2023       2.109 0.2495 Inf    1   6.312  <.0001
##  LE Year2023 / SIM Year2023       0.504 0.0447 Inf    1  -7.729  <.0001
##  LL Year2023 / NAT Year2023       3.585 0.3958 Inf    1  11.563  <.0001
##  LL Year2023 / SIM Year2023       0.857 0.0667 Inf    1  -1.988  0.4901
##  NAT Year2023 / SIM Year2023      0.239 0.0264 Inf    1 -12.981  <.0001
## 
## P value adjustment: tukey method for comparing a family of 8 estimates 
## Tests are performed on the log odds ratio scale

TSC - Figures

Early season seeded cover

ESC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year, 
##     family = binomial, data = Total_Cover_Early)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.8127  -1.5620  -0.8626   0.7913   5.1753  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.9526     0.1056 -27.955  < 2e-16 ***
## TreatmentLL   -3.3561     0.3411  -9.838  < 2e-16 ***
## TreatmentNAT  -2.0984     0.1995 -10.519  < 2e-16 ***
## TreatmentSIM  -0.3623     0.1070  -3.387 0.000707 ***
## Year2023       1.0821     0.1106   9.780  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 734.86  on 47  degrees of freedom
## Residual deviance: 247.79  on 43  degrees of freedom
## AIC: 381.05
## 
## Number of Fisher Scoring iterations: 6
## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##           LR Chisq Df Pr(>Chisq)    
## Treatment   373.50  3  < 2.2e-16 ***
## Year        106.69  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL      28.6759 9.7820 Inf    1   9.838  <.0001
##  LE / NAT      8.1535 1.6265 Inf    1  10.519  <.0001
##  LE / SIM      1.4366 0.1537 Inf    1   3.387  0.0039
##  LL / NAT      0.2843 0.1089 Inf    1  -3.284  0.0056
##  LL / SIM      0.0501 0.0172 Inf    1  -8.706  <.0001
##  NAT / SIM     0.1762 0.0360 Inf    1  -8.502  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

ESC - Figures

Fall dispersing species only

LSC - Model fitting

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
##    Data: Total_Cover_Late
## 
##      AIC      BIC   logLik deviance df.resid 
##    605.3    616.5   -296.6    593.3       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3147 -1.7714 -0.6057  1.2574  6.5425 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block  (Intercept) 0.008129 0.09016 
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.74863    0.09364 -29.352  < 2e-16 ***
## TreatmentLL   1.59358    0.09078  17.554  < 2e-16 ***
## TreatmentNAT  0.16919    0.11161   1.516  0.12953    
## TreatmentSIM  1.55200    0.09191  16.885  < 2e-16 ***
## Year2023      0.13908    0.05056   2.751  0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796                     
## TreatmntNAT -0.643  0.664              
## TreatmntSIM -0.789  0.807  0.656       
## Year2023    -0.279  0.018 -0.002  0.024
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                Chisq Df Pr(>Chisq)    
## (Intercept) 861.5598  1    < 2e-16 ***
## Treatment   566.5093  3    < 2e-16 ***
## Year          7.5683  1    0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL        0.203 0.0184 Inf    1 -17.554  <.0001
##  LE / NAT       0.844 0.0942 Inf    1  -1.516  0.4278
##  LE / SIM       0.212 0.0195 Inf    1 -16.885  <.0001
##  LL / NAT       4.155 0.3536 Inf    1  16.738  <.0001
##  LL / SIM       1.042 0.0592 Inf    1   0.732  0.8841
##  NAT / SIM      0.251 0.0216 Inf    1 -16.031  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

LSC - Figures

Plot Panel

Community Analysis

RDA Ordination

Year: 2022

Year: 2023

### fit a distance-based RDA 
#### Allows you to do a constrained ordination using non-Euclidean distance measures such as Bray-curtis dissimilarity (which is fitting for our data on plant communities)

#### it works by: 1) making the distance matrix based on your desired measure, 2) doing a PCoA, 3) taking the PCoA eigenvectors and putting them into an RDA.


### Use the capscale function (not the dbRDA since it doesn't provide species scores)
dbRDA <- capscale(inner.max.2023 ~ Treatment+Middle, data = env.data.2023,  dist="bray")

## dbRDA results
dbRDA
## Call: capscale(formula = inner.max.2023 ~ Treatment + Middle, data =
## env.data.2023, distance = "bray")
## 
##                Inertia Proportion Rank
## Total          8.39021    1.00000     
## Constrained    2.96793    0.35374    5
## Unconstrained  5.65598    0.67412   24
## Imaginary     -0.23369   -0.02785    5
## Inertia is squared Bray distance 
## Species scores projected from 'inner.max.2023' 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3   CAP4   CAP5 
## 1.3067 0.6666 0.4553 0.3562 0.1831 
## 
## Eigenvalues for unconstrained axes:
##   MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8 
## 0.7930 0.6676 0.5657 0.4643 0.4251 0.3832 0.3259 0.3169 
## (Showing 8 of 24 unconstrained eigenvalues)
## PERMANOVA results
anova(dbRDA, step = 1000, by = "terms", permu = 200)
## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = inner.max.2023 ~ Treatment + Middle, data = env.data.2023, distance = "bray")
##           Df SumOfSqs      F Pr(>F)    
## Treatment  4   2.4677 2.6178  0.001 ***
## Middle     1   0.5002 2.1227  0.012 *  
## Residual  24   5.6560                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot2 <- ordiplot(dbRDA, choices = c(1,2))

sites.long2 <- sites.long(plot2, env.data = env.data.2023)


head(sites.long2)
##   Block Treatment Year Bare_cover  Above Ground Middle       axis1      axis2
## 1     1        LE 2023          3 1852.0  114.0 1305.5 -0.46435288 -1.4255688
## 2     1        LL 2023          3 1840.0   59.0 1363.0 -0.36081765  1.2530781
## 3     1       NAT 2023          2 1854.5  127.0 1335.0 -0.04903448 -0.3763299
## 4     1       NON 2023          4 1806.0   43.5  671.0  1.04202886  0.8878276
## 5     1       SIM 2023         13 1819.5   87.5 1486.5 -1.18095045 -0.3935962
## 6     2        LE 2023         18 1786.0  196.5 1252.0  0.19736143 -1.4124790
##   labels
## 1      1
## 2      2
## 3      3
## 4      4
## 5      5
## 6      6
species.long2 <- species.long(plot2)
species.long2
##               axis1        axis2 labels
## ACHMIL -0.058008984 -0.510208190 ACHMIL
## AGRHYE -0.090320597 -0.076520219 AGRHYE
## ALLVIN -0.062351456 -0.055423685 ALLVIN
## AMBART  0.036005771 -0.127031844 AMBART
## AMBTRI  0.255047389  0.115741012 AMBTRI
## BARVUL  0.075742234  0.080996519 BARVUL
## BIDARI -0.070199953 -0.046741848 BIDARI
## BROJAP  0.087687511  0.090800725 BROJAP
## CARSPP -0.008329482 -0.025604545 CARSPP
## CERSPP -0.156534913 -0.049741117 CERSPP
## CHAFAS -0.214205742  0.070337697 CHAFAS
## CONCAN  0.254216008  0.144906241 CONCAN
## CORLAN -0.103646306 -0.238880443 CORLAN
## CORTRI -0.372322690  0.038493999 CORTRI
## CYNDAC  0.105773316  0.050351793 CYNDAC
## DESSPP  0.046739204  0.009500142 DESSPP
## ECHCRU -0.083299566 -0.019166459 ECHCRU
## ERASPE -0.138689185  0.065048505 ERASPE
## ERIANN -0.029273463  0.071860807 ERIANN
## ERYYUC -0.318825563  0.109071207 ERYYUC
## FESARU -0.041648471 -0.062619061 FESARU
## GALAPA  0.044397969 -0.108421281 GALAPA
## GALPED -0.099754338 -0.029167951 GALPED
## GERCAR  0.124078409  0.132483643 GERCAR
## GEUCAN  0.015747167 -0.009337348 GEUCAN
## HELAUT  0.056640393  0.025594942 HELAUT
## HELMOL -0.353685656  0.164837885 HELMOL
## HORPUS  0.195494739  0.220654291 HORPUS
## IPOHED  0.105773316  0.050351793 IPOHED
## IPOLAC -0.014429565  0.060524586 IPOLAC
## IVAANN -0.154551424 -0.046742336 IVAANN
## JUNINT -0.108982302  0.070549160 JUNINT
## KOEMAC -0.008746889 -0.275074186 KOEMAC
## LEPVIR  0.167659097 -0.206725676 LEPVIR
## LESCAP -0.216729478  0.145955400 LESCAP
## LIAPYC -0.121062510 -0.042119423 LIAPYC
## LONMAC -0.127721583 -0.181574047 LONMAC
## LOTCOR -0.179556651 -0.203237739 LOTCOR
## MELSPP  0.171183400  0.083139641 MELSPP
## MONFIS -0.084054969 -0.191315664 MONFIS
## MYOVER  0.306776454 -0.235079781 MYOVER
## OXADIL  0.177099525 -0.025248593 OXADIL
## PANDIC  0.018092319 -0.046244323 PANDIC
## PENDIG  0.157182736 -0.064837910 PENDIG
## PHAARU -0.016133670 -0.026788672 PHAARU
## PLALAN -0.016041685 -0.003680926 PLALAN
## PLAMAJ  0.105361035  0.004965493 PLAMAJ
## PLAVIR  0.176181476 -0.042260049 PLAVIR
## POAPRA -0.155178506 -0.044980390 POAPRA
## PYCTEN  0.033802919  0.001637247 PYCTEN
## RANABO -0.057204692  0.004765482 RANABO
## RATPIN -0.519755867  0.358576742 RATPIN
## RUDHIR -0.158255667  0.026219980 RUDHIR
## RUMCRI  0.044999745 -0.073745584 RUMCRI
## SETFAB  0.319110355  0.170382494 SETFAB
## SETPUM  0.565939547  0.007443133 SETPUM
## SIDSPI  0.102095417  0.041813388 SIDSPI
## SILANT  0.106599037  0.050853680 SILANT
## SOLCAN  0.083506391  0.036817574 SOLCAN
## SOLCAR  0.154563253  0.069255882 SOLCAR
## SORNUT -0.458619292  0.200809834 SORNUT
## SPHOBT -0.127212271 -0.014392910 SPHOBT
## SPOHET -0.121062510 -0.042119423 SPOHET
## SYMPIL  0.078367533 -0.137094939 SYMPIL
## TAROFF  0.119137849 -0.062499146 TAROFF
## TRIPRA -0.121062510 -0.042119423 TRIPRA
## TRIREP  0.176311277 -0.181646697 TRIREP
## VERARV  0.412681869  0.226543389 VERARV
## VERURT -0.108982302  0.070549160 VERURT
## VIOPED  0.044027264 -0.176350164 VIOPED
axis.long2 <- axis.long(dbRDA, choices = c(1,2))
axis.long2 
##   axis     ggplot        label
## 1    1 xlab.label CAP1 (15.6%)
## 2    2 ylab.label  CAP2 (7.9%)
BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 12, family="Arial"),
axis.text = element_text(size = 10, colour = "gray25"),
axis.title = element_text(size = 14, colour = "gray25"),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())
plotgg2 <- ggplot() + 
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axis.long2[1, "label"]) +
    ylab(axis.long2[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_point(data=sites.long2, 
               aes(x=axis1, y=axis2, colour=Treatment, shape=Treatment), 
               size=5) +
    geom_point(data=species.long2, 
               aes(x=axis1, y=axis2)) +
    BioR.theme +
    ggsci::scale_colour_npg() +
    coord_fixed(ratio=1)
spec.envfit <- envfit(plot2, env=inner.max.2023)
spec.data.envfit <- data.frame(r=spec.envfit$vectors$r, p=spec.envfit$vectors$pvals)
species.long2 <- species.long(plot2, spec.data=spec.data.envfit)
species.long2
##                   r     p        axis1        axis2 labels
## ACHMIL 0.5072454500 0.001 -0.058008984 -0.510208190 ACHMIL
## AGRHYE 0.0639828837 0.407 -0.090320597 -0.076520219 AGRHYE
## ALLVIN 0.0377623714 0.611 -0.062351456 -0.055423685 ALLVIN
## AMBART 0.0599078832 0.441  0.036005771 -0.127031844 AMBART
## AMBTRI 0.2392522931 0.028  0.255047389  0.115741012 AMBTRI
## BARVUL 0.1303343327 0.143  0.075742234  0.080996519 BARVUL
## BIDARI 0.0168963769 0.820 -0.070199953 -0.046741848 BIDARI
## BROJAP 0.0651875300 0.398  0.087687511  0.090800725 BROJAP
## CARSPP 0.0248129138 0.744 -0.008329482 -0.025604545 CARSPP
## CERSPP 0.0287982535 0.672 -0.156534913 -0.049741117 CERSPP
## CHAFAS 0.2000449585 0.033 -0.214205742  0.070337697 CHAFAS
## CONCAN 0.2824451407 0.010  0.254216008  0.144906241 CONCAN
## CORLAN 0.2291601219 0.022 -0.103646306 -0.238880443 CORLAN
## CORTRI 0.3719211480 0.001 -0.372322690  0.038493999 CORTRI
## CYNDAC 0.1014468817 0.309  0.105773316  0.050351793 CYNDAC
## DESSPP 0.0150879971 0.912  0.046739204  0.009500142 DESSPP
## ECHCRU 0.0465040615 0.694 -0.083299566 -0.019166459 ECHCRU
## ERASPE 0.1031003723 0.251 -0.138689185  0.065048505 ERASPE
## ERIANN 0.0522202337 0.531 -0.029273463  0.071860807 ERIANN
## ERYYUC 0.2378151015 0.023 -0.318825563  0.109071207 ERYYUC
## FESARU 0.0008621138 0.991 -0.041648471 -0.062619061 FESARU
## GALAPA 0.0498040088 0.692  0.044397969 -0.108421281 GALAPA
## GALPED 0.0297373414 0.647 -0.099754338 -0.029167951 GALPED
## GERCAR 0.1152105947 0.187  0.124078409  0.132483643 GERCAR
## GEUCAN 0.0172861958 0.853  0.015747167 -0.009337348 GEUCAN
## HELAUT 0.0128013368 0.917  0.056640393  0.025594942 HELAUT
## HELMOL 0.4037999010 0.002 -0.353685656  0.164837885 HELMOL
## HORPUS 0.2497223863 0.011  0.195494739  0.220654291 HORPUS
## IPOHED 0.1014468817 0.309  0.105773316  0.050351793 IPOHED
## IPOLAC 0.0269072643 0.731 -0.014429565  0.060524586 IPOLAC
## IVAANN 0.1299775945 0.119 -0.154551424 -0.046742336 IVAANN
## JUNINT 0.0994951469 0.326 -0.108982302  0.070549160 JUNINT
## KOEMAC 0.2098574477 0.025 -0.008746889 -0.275074186 KOEMAC
## LEPVIR 0.1922656027 0.075  0.167659097 -0.206725676 LEPVIR
## LESCAP 0.2672247308 0.003 -0.216729478  0.145955400 LESCAP
## LIAPYC 0.1084649946 0.240 -0.121062510 -0.042119423 LIAPYC
## LONMAC 0.2098758530 0.038 -0.127721583 -0.181574047 LONMAC
## LOTCOR 0.2377141991 0.017 -0.179556651 -0.203237739 LOTCOR
## MELSPP 0.1105157989 0.198  0.171183400  0.083139641 MELSPP
## MONFIS 0.1949009232 0.044 -0.084054969 -0.191315664 MONFIS
## MYOVER 0.2442572475 0.024  0.306776454 -0.235079781 MYOVER
## OXADIL 0.1444608064 0.091  0.177099525 -0.025248593 OXADIL
## PANDIC 0.0021488040 0.981  0.018092319 -0.046244323 PANDIC
## PENDIG 0.0933758114 0.281  0.157182736 -0.064837910 PENDIG
## PHAARU 0.0089886654 0.966 -0.016133670 -0.026788672 PHAARU
## PLALAN 0.0001244903 0.999 -0.016041685 -0.003680926 PLALAN
## PLAMAJ 0.0481384833 0.624  0.105361035  0.004965493 PLAMAJ
## PLAVIR 0.1187269285 0.158  0.176181476 -0.042260049 PLAVIR
## POAPRA 0.1689308407 0.060 -0.155178506 -0.044980390 POAPRA
## PYCTEN 0.0148643439 0.936  0.033802919  0.001637247 PYCTEN
## RANABO 0.0176453416 0.799 -0.057204692  0.004765482 RANABO
## RATPIN 0.6577890834 0.001 -0.519755867  0.358576742 RATPIN
## RUDHIR 0.0591117238 0.429 -0.158255667  0.026219980 RUDHIR
## RUMCRI 0.0039828781 0.960  0.044999745 -0.073745584 RUMCRI
## SETFAB 0.3549839356 0.001  0.319110355  0.170382494 SETFAB
## SETPUM 0.7069156513 0.001  0.565939547  0.007443133 SETPUM
## SIDSPI 0.0813242368 0.369  0.102095417  0.041813388 SIDSPI
## SILANT 0.0406318031 0.748  0.106599037  0.050853680 SILANT
## SOLCAN 0.0370632877 0.769  0.083506391  0.036817574 SOLCAN
## SOLCAR 0.1279225355 0.124  0.154563253  0.069255882 SOLCAR
## SORNUT 0.5060084303 0.001 -0.458619292  0.200809834 SORNUT
## SPHOBT 0.0259085617 0.751 -0.127212271 -0.014392910 SPHOBT
## SPOHET 0.1084649946 0.240 -0.121062510 -0.042119423 SPOHET
## SYMPIL 0.1437469470 0.123  0.078367533 -0.137094939 SYMPIL
## TAROFF 0.0595860076 0.427  0.119137849 -0.062499146 TAROFF
## TRIPRA 0.1084649946 0.240 -0.121062510 -0.042119423 TRIPRA
## TRIREP 0.2735959378 0.015  0.176311277 -0.181646697 TRIREP
## VERARV 0.3327591785 0.007  0.412681869  0.226543389 VERARV
## VERURT 0.0994951469 0.326 -0.108982302  0.070549160 VERURT
## VIOPED 0.2171418847 0.002  0.044027264 -0.176350164 VIOPED
species.long3 <- species.long2[species.long2$p < 0.05, ]
species.long3
##                r     p        axis1        axis2 labels
## ACHMIL 0.5072454 0.001 -0.058008984 -0.510208190 ACHMIL
## AMBTRI 0.2392523 0.028  0.255047389  0.115741012 AMBTRI
## CHAFAS 0.2000450 0.033 -0.214205742  0.070337697 CHAFAS
## CONCAN 0.2824451 0.010  0.254216008  0.144906241 CONCAN
## CORLAN 0.2291601 0.022 -0.103646306 -0.238880443 CORLAN
## CORTRI 0.3719211 0.001 -0.372322690  0.038493999 CORTRI
## ERYYUC 0.2378151 0.023 -0.318825563  0.109071207 ERYYUC
## HELMOL 0.4037999 0.002 -0.353685656  0.164837885 HELMOL
## HORPUS 0.2497224 0.011  0.195494739  0.220654291 HORPUS
## KOEMAC 0.2098574 0.025 -0.008746889 -0.275074186 KOEMAC
## LESCAP 0.2672247 0.003 -0.216729478  0.145955400 LESCAP
## LONMAC 0.2098759 0.038 -0.127721583 -0.181574047 LONMAC
## LOTCOR 0.2377142 0.017 -0.179556651 -0.203237739 LOTCOR
## MONFIS 0.1949009 0.044 -0.084054969 -0.191315664 MONFIS
## MYOVER 0.2442572 0.024  0.306776454 -0.235079781 MYOVER
## RATPIN 0.6577891 0.001 -0.519755867  0.358576742 RATPIN
## SETFAB 0.3549839 0.001  0.319110355  0.170382494 SETFAB
## SETPUM 0.7069157 0.001  0.565939547  0.007443133 SETPUM
## SORNUT 0.5060084 0.001 -0.458619292  0.200809834 SORNUT
## TRIREP 0.2735959 0.015  0.176311277 -0.181646697 TRIREP
## VERARV 0.3327592 0.007  0.412681869  0.226543389 VERARV
## VIOPED 0.2171419 0.002  0.044027264 -0.176350164 VIOPED
species.long3 <- species.long3 %>% 
  rename(SPP6 = "labels")

species.long3 <- left_join(species.long3, species_list_PE, by = "SPP6")
vectors.envfit <- envfit(plot2, env=env.data.2023)
vectors.long3 <- vectorfit.long(vectors.envfit)
vectors.long3
##                vector      axis1      axis2          r     p
## Block           Block  0.8032309  0.5956678 0.02300445 0.729
## Bare_cover Bare_cover  0.9653136 -0.2610932 0.01682397 0.806
## Above           Above -0.4497822 -0.8931382 0.01865029 0.786
## Ground         Ground -0.8836873 -0.4680778 0.06824846 0.409
## Middle         Middle -0.8443056 -0.5358620 0.25754975 0.019
ordiplot(dbRDA, choices=c(1,2))
Treatment.ellipses <- ordiellipse(plot2, groups=env.data.2023$Treatment, display="sites", kind="sd")

Treatment.ellipses.long2 <- ordiellipse.long(Treatment.ellipses, grouping.name="Treatment")

All

PERMANOVA

PERMANOVA analysis revealed that treatment did have a significant effect on plant community composition. Treatment also explained ~ 20% of the variation in the plant community. Post-hoc analysis shows that LE, NON, and NAT treatments were not significantly different from each other but did differ in plant community composition compared to the SIM and LL seeding treatments. Additionally, the SIM and LL treatments were only marginally different (p = 0.053). Based on the NMDS above, the LL plant community appears as a subset of the SIM plant community. Despite receiving the same pool of seeded species, species added to the plots later in the growing season did not establish as well. Barriers to establishment also appeared higher for spring/summer dispersing species compared to fall dispersing species, particularly when added without priority.

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = inner.max ~ Treatment + Year, data = inner_wide, permutations = 1000, method = "bray")
##           Df SumOfSqs      R2      F   Pr(>F)    
## Treatment  4   3.1762 0.17687 3.3089 0.000999 ***
## Year       1   1.8224 0.10148 7.5940 0.000999 ***
## Residual  54  12.9588 0.72164                    
## Total     59  17.9574 1.00000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## $parent_call
## [1] "inner.max ~ Treatment + Year , strata = Null , permutations 999"
## 
## $LE_vs_LL
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.9839 0.14313 4.2642  0.001 ***
## Year       1   1.0450 0.15202 4.5291  0.001 ***
## Residual  21   4.8456 0.70486                  
## Total     23   6.8745 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.4410 0.06307 1.6582  0.039 *  
## Year       1   0.9656 0.13812 3.6310  0.001 ***
## Residual  21   5.5845 0.79881                  
## Total     23   6.9911 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.6830 0.10287 2.8438  0.001 ***
## Year       1   0.9126 0.13746 3.8000  0.001 ***
## Residual  21   5.0433 0.75967                  
## Total     23   6.6389 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.7381 0.10692 2.9548  0.001 ***
## Year       1   0.9197 0.13323 3.6819  0.001 ***
## Residual  21   5.2454 0.75986                  
## Total     23   6.9031 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.7043 0.10632 2.8839  0.001 ***
## Year       1   0.7917 0.11951 3.2417  0.001 ***
## Residual  21   5.1287 0.77418                  
## Total     23   6.6247 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   1.1298 0.17501 5.2743  0.001 ***
## Year       1   0.8276 0.12820 3.8636  0.001 ***
## Residual  21   4.4985 0.69680                  
## Total     23   6.4560 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.5752 0.09413 2.6448  0.001 ***
## Year       1   0.9681 0.15842 4.4511  0.001 ***
## Residual  21   4.5672 0.74744                  
## Total     23   6.1105 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.4520 0.07021 1.8176  0.026 *  
## Year       1   0.7635 0.11859 3.0701  0.001 ***
## Residual  21   5.2222 0.81119                  
## Total     23   6.4377 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.9454 0.13241 3.6145  0.001 ***
## Year       1   0.7019 0.09831 2.6836  0.003 ** 
## Residual  21   5.4929 0.76928                  
## Total     23   7.1403 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NON_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   1.2878 0.18695 5.5308  0.001 ***
## Year       1   0.7109 0.10321 3.0533  0.003 ** 
## Residual  21   4.8897 0.70984                  
## Total     23   6.8884 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Weather data

# Precipitation



month_tot <- weather %>% 
  filter(Year != "2021") %>% 
  group_by(Month, Year) %>% 
  summarize(tot_month_PRP_mm = sum(Total_PRP_mm))
## `summarise()` has grouped output by 'Month'. You can override using the
## `.groups` argument.
year_tot <- weather %>% 
  filter(Year != "2021") %>% 
  group_by(Year) %>% 
  summarize(tot_PRP_mm = sum(Total_PRP_mm))

library(mosaic)
## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## 
## The following objects are masked from 'package:rstatix':
## 
##     cor_test, prop_test, t_test
## 
## The following object is masked from 'package:lmerTest':
## 
##     rand
## 
## The following object is masked from 'package:lme4':
## 
##     factorize
## 
## The following object is masked from 'package:Matrix':
## 
##     mean
## 
## The following object is masked from 'package:cowplot':
## 
##     theme_map
## 
## The following objects are masked from 'package:goeveg':
## 
##     deg2rad, rad2deg
## 
## The following object is masked from 'package:permute':
## 
##     shuffle
## 
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     stat
## 
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## 
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
favstats(year_tot$tot_PRP_mm)
##      min      Q1 median       Q3      max    mean       sd  n missing
##  665.988 848.995 934.72 1035.304 1148.842 927.735 143.0993 10       0
# Average air temperature


library(mosaic)

favstats(weather$Average_Air_Temp_C)
##    min       Q1   median       Q3      max     mean       sd    n missing
##  -21.5 4.277778 13.33333 21.55556 32.66667 12.39584 10.45545 4018       0